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Abstract 

Computational  methods  for  simulation  of  the  tunneling  stage  of  penetration  into  semi-infinite 
compressible  geological  targets  of  a  hard  (rigid)  penetrator  for  intermediate  impact  velocities 
(typically,  below  1000  m/s)  were  developed.  A  compressible  viscoplastic  fluid  constitutive  equation 
that  captures  the  combined  effects  of  high-strain  rate  and  high-pressure  (confinement)on  yielding 
was  developed.  To  account  for  the  experimentally  observed  characteristics  of  the  response  at 
higli-pressures  the  hypothesis  of  a  locking  medium  was  adopted  (i.e.  the  density  cannot  exceed 
a  critical  limit).  To  simulate  steady-state  penetration,  a  mixed  finite-element  and  finite- volume 
strategy  is  developed.  Specifically,  the  variational  inequality  for  the  velocity  field  is  discretized 
using  the  finite  element  method  and  a  finite  volume  method  is  adopted  for  the  density  equation.  To 
solve  the  velocity  problem  a  decomposition-coordination  formulation  coupled  with  the  augmented 
lagrangian  method  is  used.  This  approach  is  accurate  in  detecting  the  visco-plastic  regions  and 
permit  us  to  handle  the  locking  condition.  The  ability  of  the  proposed  model  to  accurately  describe 
(i)  the  density  changes  in  the  target  around  the  penetration  tunnel,  ii)  the  shape  and  location 
of  rigid/plastic  boundary,  iii)  the  directional  damage/fracture  ahead  of  the  penetrator.  Moreover 
the  fact  that  fracture  occurs  ahead  of  the  penetrator  and  along  planes  which  are  not  symmetric 
with  respect  to  the  penetrator  centerline,  can  explain  trajectory  instabilities. 
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1  Introduction 


Kinetic  energy  penetration  phenomena  are  of  interest  in  a  variety  of  applications  ranging  from 
terminal  ballistics  to  protection  of  spacecraft  due  to  meteoroid  impact,  containment  of  high 
mass  or  high  velocity  debris  due  to  accidents  or  high  rate  energy  release,  design  of  hardened 
protective  facilities,  erosion  and  fracture  of  solids  due  to  impact,  etc.  (Zukas  [32]).  The 
occurrence  of  multiple  phenomena  in  the  target  such  as  localization,  plasticity,  anisotropic 
damage,  fragmentation  pushes  the  limits  of  existing  modelling  and  computational  capabilities 
for  description  of  the  target  response. 

Since  the  deformation  rates  are  very  large  (of  the  order  of  105s-1),  a  fluid  formulation  in 
eulerian  description  could  describe  the  main  features  of  the  target  deformation.  Because  shear 
flow  takes  place  only  if  a  certain  threshold  is  surpassed,  a  non-newtonian  fluid  model  is  not 
appropriate.  Since  yielding  cannot  be  neglected,  viscoplastic  fluid  models  such  as  Bingham 
have  been  used  [2,  9,  17,  18]  in  studies  concerning  penetration  into  metalic  targets.  It  is  to 
be  noted  that  such  models  are  suitable  for  the  description  of  the  high-strain  rate  behavior 
of  metals  for  which  yielding  is  insensitive  to  pressure  and  consequently  are  incompressible. 
As  opposed  to  metals,  in  geological  materials  yielding  is  pressure-dependent  (e.g.  [27]).  For 
these  materials,  compressibility  is  not  a  second  order  property  as  for  non-newtonian  fluids. 
Indeed,  the  geologic  and  cementitious  materials  have  a  compaction  rate  up  to  30%  and  their 
mechanical  properties  (viscosity,  yield  limit,  etc)  depend  significantly  on  density  (4  to  400 
times  increase  in  yield  limit  for  quasi-static  and  dynamic  conditions,  respectively).  Thus, 
compressible  viscoplastic  fluid  models  should  be  considered. 

Only  a  limited  number  of  2-D  and  3-D  calculations  of  penetration  into  cementitious  or 
geologic  materials  for  low  to  intermediate  velocity  impacts  have  been  reported  in  the  open  lit¬ 
erature.  Finite-element  and  finite-difference  models  employing  both  Lagrangian  and  Eulerian 
formulations  have  been  used.  Each  of  these  methods  has  advantages  and  disadvantages  as 
a  tool  to  model  penetration.  For  example,  with  the  finite-element  method  using  Lagrangian 
formulation,  one  needs  to  perform  a  remeshing  whenever  the  mesh  becomes  highly  distorted 
as  the  projectile  advances  in  the  target.  A  finite-difference  method  generally  utilizes  an  Eule¬ 
rian  description  in  which  the  mesh  is  fixed  in  space,  and  thus  does  not  distort.  The  moving 
boundary  between  the  penetrator  and  the  surrounding  target  is  not  sharp,  and  special  inter¬ 
face  tracking  methods  are  required  (e.g.  level  set  methods).  Mesh-less  methods  have  also  been 
proposed  for  tracking  moving  fronts  of  discontinuities. 

In  discrete  element  models  (DEM)  the  material  is  treated  as  heterogeneous  at  the  macro¬ 
scale.  Specifically,  the  material  is  discretized  into  individual  elements  or  blocks,  which  are 
allowed  to  interact  as  stress  is  applied  [23].  Applications  of  DEM  to  modelling  the  dynamic 
behavior  of  concrete  have  been  reported  [6,  10]  and  a  very  good  agreement  with  experiments 
was  obtained.  However,  DEM  is  not  applicable  to  large  scale  problems  because  the  computa¬ 
tion  time  required  to  solve  even  simple,  routine  problems  can  be  excessive. 

It  is  to  be  noted  that  the  material  models  used  in  conjunction  with  the  aforementioned 
methods  exhibit  limited  features  such  as  pressure-dependent  yield  surfaces  (e.g.  Schwer  et  al. 
[27]),  strain- rate  dependent  yield  surfaces  (e.g.  Batra  [2]),  pressure  and  strain-rate  dependent 
yield  surfaces  (e.g.  Adarnik  and  Matejovic,  [1]),  simple  equations  of  state  for  porous  materials 
models  (e.g.  Tipton  [31]).  Most  of  the  material  models  are  incompressible. 

A  new  model  for  describing  steady-state  penetration  in  geological  or  geologically  derived 
materials  is  presented  in  this  paper  which  is  outlined  as  follows.  To  account  for  the  combined 
effects  of  high-pressure  and  high  strain-rate  on  the  flow  behavior,  a  new  constitutive  equation 
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is  proposed  in  section  2.  The  focus  is  on  capturing  the  increase  in  the  yielding  limit  with 
the  degree  of  compaction  as  well  as  the  observed  decrease  in  strain  rate  sensitivity.  Since 
in  the  range  of  impact  velocities  of  interest,  the  target  material  displays  both  solid-like  and 
fluid-like  behavior,  the  material  model  will  be  obtained  by  superposing  a  rigid-plastic  solid  to 
a  compressible  viscous  fluid. 

The  problem  statement  of  steady  state  penetration  and  its  variational  formulations  are 
presented  in  section  3.  A  mixed  finite-element  and  finite- volume  strategy  is  developed  in 
section  4.  Specifically,  the  variational  inequality  for  the  velocity  field  is  discretized  using  the 
finite  element  method  and  a  finite  volume  method  is  adopted  for  the  density  equation.  To  solve 
the  velocity  problem  a  mixed  formulation  with  the  augmented  lagrangian  method  is  used.  The 
model  is  further  applied  to  axisymmetric  penetration  into  concrete  in  section  5.  The  material 
parameters  are  found  following  the  recent  experimental  data  obtained  by  Schmidt  [29].  In 
order  to  capture  sharply  the  shape  of  the  visco-plastic  zone,  we  have  used  an  anisotropic  mesh 
generator.  The  model  predicts  that  around  the  penetrator,  a  fully  compacted  state  is  achieved, 
the  maximum  compaction  being  in  the  nose  zone.  Fracture  occurs  ahead  of  the  penetrator 
and  along  planes  which  can  be  not  symmetric  with  respect  to  the  penetrator  centerline.  This 
induced  anisotropy  can  explain  trajectory  deviations 


2  Material  modeling 

In  the  hypervelocity  impact  range,  the  flow  behavior  of  the  target  could  be  described  accurately 
by  classic  fluid  type  constitutive  equations.  At  low-to- intermediate  impact  velocity  (below  1000 
m/s),  the  impacted  medium  may  display  both  solid  and  fluid-like  properties.  Thus,  model  for 
fluids  with  yield  limit  are  generally  used.  Since  yielding  of  metals  is  insensitive  to  pressure, 
these  models  are  incompressible  (e.g.  Bingham  type  constitutive  equation)  [2,  9,  17,  18]).  In 
contrast,  for  geological  materials  yielding  is  pressure-dependent  (e.g.  [27]).  The  effect  of  the 
density  on  flow  cannot  be  neglected,  since  compressibility  is  not  a  second-  order  property  as 
for  non-newtonian  fluids.  Thus,  a  compressible  Bingham  type  model  with  yielding  dependent 
on  the  actual  density  (or  compaction  level)  needs  to  be  considered. 

Let  denote  by  u  the  material  velocity,  by  D  the  rate  of  deformation  tensor  and  by  D'  = 
D  —  |(trZ))/3  its  deviator 

D  =  D(u)  =  -(Vu  +  VTu),  D'  =  D'[u)  =  D(u)  -  -divu/3.  (1) 

2  3 

Also,  the  Cauchy  stress  tensor  is  denoted  by  er  and  its  deviator  by  <j'  =  a  —  tr  er/3.  In  contrast 
to  a  Navier-Stokes  fluid,  a  classic  Bingham  fluid  (see  [3,  22,  8]),  can  sustain  a  shear  stress 
even  at  rest  and  it  starts  to  flow  only  if  the  applied  forces  exceed  an  yield  limit  k.  To  account 
for  the  effect  of  compaction  on  the  deviatoric  response  of  cementitious  or  geologic  materials, 
in  [4]  the  yield  limit  k  was  considered  to  be  a  function  of  the  current  density.  In  this  paper, 
we  consider  an  extension  (see  [11])  of  the  Bingham  model  which  is  obtained  by  superposing  a 
rigid-plastic  solid  to  a  compressible  viscous  fluid,  i.e. 

a  =  S  +  [-p(p)  +  X{p){tT  D)\I3  +  2rj{p)D 

F(S,p)<  0 
D'  =  pdsT 
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(2) 

(3) 

(4) 


where  the  tensor  S  is  the  part  of  the  stress  which  describes  the  plastic  properties  of  the 
material,  p  is  the  density  while  77,  A  >  0,  are  viscosity  coefficients ,  T  is  the  yield  function  and 
p  is  a  scalar  function  such  that 

ji(t)  =  0  if  tF(S,p)  <  0  or  tF(S,p)  =  0  and  dstF(S,p)  :  S  <  0 

ffit)  >0  if  F(S,  p)  =  0  and  %F(S,  p)  :  S  =  0  '  " ' 


Note  that  the  viscosity  coefficients  77,  A  as  well  as  the  yield  function  T  are  considered  to  be 
functions  of  the  current  density  p  or  alternatively  on  the  compaction  factor  c  =  p/po  —  1- 
Further,  it  will  be  assumed  that  the  Mises  condition  is  satisfied,  i.e. 

F(S)  =  \S'\2-K2(p),  (6) 

where  n2(p)  is  the  yield  limit  in  shear  which  is  considered  to  be  a  function  of  the  current  density. 
Experimental  studies  of  the  dynamic  behavior  of  cementitious  materials  have  indicated  a  strong 
dependence  of  yielding  and  subsequent  flow  on  density  (see  for  instance  experimental  data  on 
cementitious  materials  reported  in  [21,  26,  28]). 

The  material  function  p(p),  describes  the  volumetric  response  of  the  material.  Under  hy¬ 
drostatic  conditions,  most  cementitious  materials  show  a  highly  non-linear  pressure- volumetric 
strain  response,  the  reversible  decrease  in  volume  being  very  small.  The  experimental  observa¬ 
tions  also  suggest  that  in  the  high-pressure  regime,  a  very  large  increase  in  pressure  is  necessary 
in  order  to  produce  even  a  very  small  change  in  density.  Thus,  the  hypothesis  of  a  ” locking 
medium”  can  be  made,  i.e.  the  density  cannot  exceed  a  critical  value.  This  critical  density 
p*,  called  locking  density,  corresponds  to  a  state  in  the  material  when  all  the  pores  and  cracks 
are  closed.  The  pressure  level  at  which  this  density  is  first  reached,  called  locking  pressure,  is 
denoted  by  p*  =  p{p*)-  Hence, 

(  P  =  P(p),  if  P  <  P* 

(7) 

{  P  >P*,  if  P  =  P* 


Since  during  unloading  the  reversible  decrease  of  volume  is  very  small,  it  can  be  neglected 
and  therefore  the  unloading  process  is  rigid.  Hence, 


tr  D  <  0, 


tr  S  =  0,  if  trD  <  0 

tr  S  >  0,  if  tr  D  =  0. 


(8) 


Following  the  procedure  used  to  construct  the  classical  Bingham  fluid  (see  [8] ) ,  from  equations 
(2)  to  (5)  we  deduce  the  relationship  between  the  deviatoric  part  of  the  rate  of  deformation 
and  the  stress  deviator: 


D'  = 


1 

2  ffip)  +  a 


if  \a'\  >  n(p) 
if  | cr'  |  <  k(p). 


(9) 


Next, 

by 


by  inverting  the  constitutive  equation  (9)  we  obtain  that  the  shear  response  is  governed 


a' =  2r]{p)D' +  k{p)^~  if  \D'\  +  0 
| <r' |  <  k(p)  if  \D'\  =  0. 


(10) 
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The  formulation  of  the  model  is  completed  by  providing  the  equations  governing  the  response 
under  hydrostatic  conditions.  From  (2)  and  (8) 


trcr  =  —3 p(p)  +  (3A (p)  +  2rj(p))tr  D  if  trD  <  0 
trcr  >  —3 p(p)  if  tr  D  =  0. 


(11) 


Subsequently  we  will  refer  to  (10)-(11)  as  the  constitutive  equations  of  the  compressible 
Binghanr-type  material.  Note  that  the  classical  Bingham  fluid  (1)  is  recovered  if  in  (10)  the 
incompressibility  condition  (i.e.  diva  =  0)  is  imposed. 


3  Statement  of  the  problem 

In  this  section  we  present  the  equations  governing  the  steady-state  motion  of  concrete,  de¬ 
scribed  by  the  constitutive  equation  (10)-(11),  over  a  rigid  penetrator  fully  embedded  in  the 
target,  represented  by  a  domain  Pci3  with  a  smooth  boundary  &D. 


E 


Figure  1:  A  schematic  representation  of  the  domain  V. 


The  domain  T>  is  assumed  to  be  the  whole  space  M3  without  the  penetrator  V  and  the 
infinite  tunnel  T  behind  it.  8qD  where  the  velocities  V  are  prescribed  will  be  the  infinity  of 
the  domain  V;  d\D  is  the  boundary  of  the  tunnel  while  82V,  is  the  part  of  the  boundary  where 
there  is  frictional  contact  with  the  projectile.  However,  post-test  observations  indicate  that 
the  tunnel  is  of  the  order  of  the  projectile  diameter.  For  both  grout  and  concrete  targets  a 
change  in  density  in  the  region  around  the  penetration  tunnel  radially  outward  from  the  edge 
of  the  tunnel  to  a  distance  of  1  —  1.5  projectile  diameters  was  reported  (see  Jones  et  al  [19]). 
Hence,  the  domain  affected  by  the  impact  event  is  bounded  and  thus  V  will  be  restricted 
to  the  domain  shown  Figure  1.  To  enable  reasonable  computational  effort  and  still  ensure 
that  the  boundary  conditions  at  infinity  are  accurately  described,  we  limit  the  extent  of  this 
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domain  to  5  projectile  radii.  The  projectile  is  considered  to  be  rigid  and  axisymmetric.  Since, 
at  striking  velocities  up  to  1200  m/s,  penetration  paths  are  relatively  straight  and  stable  with 
regard  to  the  original  shotline,  the  problem  could  be  considered  axisymmetric  with  respect  to 
the  projectile  centerline  0 z. 

The  momentum  balance  law  in  the  Eulerian  coordinates  reads 


p(u  ■  V)w  —  div  a  =  pf  in  V, 


(12) 


where  p  =  p(x)  >  po  >  0  is  the  mass  density  distribution  and  /  denotes  the  body  forces.  The 
continuity  equation  is 

di v(pu)  =0  in  V.  (13) 

To  the  system  of  equations  (10)-(ll)and  (12)we  associate  the  following  boundary  conditions: 


u  =  V  on  OqD,  an  =  0  on  d\  V,  (14) 

where  n  stands  for  the  outward  unit  normal  on  dT>,  un  =  u  ■  n  is  normal  velocity,  at  = 
a  —  (a  ■  n)n  stands  for  the  tangential  stress  and  V  is  the  imposed  velocity. 

A  slip-dependent  frictional  contact  is  assumed  on  the  boundary  d-fiD 


Un  —  0) 


| O’* |  <  nWn], 

l  l  Ut 

&t  —  0"n  |  |  i 

ml 


(15) 


where  an  =  an  •  n  is  the  normal  stress,  Ut  =  u  —  unn  is  the  tangential  velocity  and  p  is  the 
friction  coefficient.  According  to  (15),  the  the  tangential  (friction)  stress  is  bounded  by  the 
normal  stress  multiplied  by  the  friction  coefficient  p.  If  such  a  limit  is  not  attained  sliding 
cannot  occur;  otherwise  the  friction  stress  is  opposite  to  the  slip  rate.  The  friction  coefficient 
will  be  considered  variable  during  the  slip.  Experimental  observations  indicate  that  the  friction 
coefficient  depends  on  the  slip  rate  |wj|  i.e.  p  =  p(\ut\).  The  simplest  law  of  variation  of  p  on 
the  slip  rate  is  a  discontinuous  jump  from  a  “static”  value  (for  \ut]  =  0)  down  to  a  “dynamic” 
or  “kinetic”  value  (for  \ut\  0).  In  this  work,  we  will  consider  a  smooth  and  decreasing 
function  p  =  p(\ut\)  of  the  slip  rate. 

The  boundary  condition  for  the  conservation  of  mass  equation  is: 


P  =  po  on  d PV  (16) 

i.e.  a  given  density  is  prescribed  on  dpT>  C  d T>. 

Setting 

k={r£  i/1(T))3;u  =  V  on  doV,  vn  =  0  on  V  =  {v  e  V;  div  v  <  0  in  V  } 

the  kinematic  admissible  set,  the  variational  formulation  for  the  velocity  field  u  e  V  is: 

f  p(u  ■  V)u  ■  (v  —  u)  +  f  2r](p)D'(u)  :  (. D'(v )  —  D\u))+ 

Jv  Jv 

[  k{p){\D\v)\-\D'{u)\)+  f  [(A(p)  +  !??(p))divw  -p(p)]div(«  -  u) 

Jv  Jv  6 

+  [  MM)KI(M  -  |«t|)  >  [  pf-{v-u).  (17) 

Jd2V  Jv 

for  all  v  £  V. 

Thus,  the  problem  of  the  flow  of  compressible  visco- plastic  Binghanr-type  fluid  becomes: 
Find  the  velocity  field  u  and  the  mass  density  field  p  such  that  equations  (13),  and  (17) 
hold. 
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4  Numerical  approach 

In  this  section,  we  present  the  numerical  approach  adopted  for  solving  the  problem  of  the 
stationary  flow  of  a  rigid-viscoplastic  Bingham  type  material. 

The  algorithm  consists  of  solving  alternatively  the  variational  inequality  (17)  for  the  ve¬ 
locity  field  and  the  continuity  equation  (13)  for  the  density  field.  More  precisely,  we  shall 
distinguish  two  problems  :  the  ” velocity  problem”  and  the  the  ” density  problem”.  For  the 
velocity  problem  we  assume  that  p  and  the  distributions  of  p.  77,  A  and  k  are  given,  and  we 
find  u  6  V,  the  solution  of  (17).  The  density  problem  consists  in  finding  the  density  field  p 
solution  of  (13)  assuming  that  u  is  given. 

As  it  follows  from  the  next  subsection  the  velocity  problem  will  be  solved  by  an  iterative 
procedure.  At  each  iteration  k  we  get  a  velocity  field  uk  which  will  be  used  to  solve  the  density 
problem 

di v{pkuk)  =0  in  V,  pk  =  p0  on  dpV  (18) 

and  getting  pk  to  update  the  density  field. 


4.1  The  velocity  problem 

Suppose  that  V  is  discretized  by  using  a  family  of  triangulations  (7 h)h  made  of  finite  elements 
of  degree  2  where  h  >  0  is  the  discretization  parameter  representing  the  greatest  diameter  of 
a  triangle  in  Tk.  The  finite  element  space  V/,.,  which  is  an  internal  appoximation  of  V  reads: 

Vh  =  [vh-  vh  G  ( C(V ))3,  vh\T  g  (P2(T))3  VT  g  Th ,  vh  =  V  on  d0V ,  vhn  =  0  on  d2V  }, 

where  C(V)  stands  for  the  space  of  continuous  functions  on  V  and  Pi(T )  represents  the  space 
of  polynomial  functions  of  degree  i  on  T. 

Let  Vh  =  V  f)  Vh ■  Then,  the  velocity  problem  is  discretized  by  considering  Uh  G  Vh  which 
satisfies  (17)  for  all  v  G  Vh-  In  order  to  simplify  the  notations  we  shall  omit  in  all  this 
subsection  the  indexes  h. 

For  all  w  eVh  let  Jw  :  V),,  — >  M  be  given  by 

Jwiy)  =  [  v\D'(v)\2  +  f  (div  v)2  +  f  n\D'(v)\  +  f  p,(\wt\)\an(w)\\vt\ 

Jv  Jv  \z  6  J  Jv  Jd2T> 


+  /  p(w-'V)w-v—  j  pdivv—  j  pf-v. 

Jv  . 

Hence,  the  velocity  problem  can  be  written  as 


lv 


lv 


U  G  Vh,  ju(u)  =  inf  Ju(v) 
vevh 


(19) 


(20) 


The  following  iterative  algorithm  is  proposed  :  given  uk  1  G  Vh  find  uk  G  Vh ,  the  solution  of 
the  following  minimization  problem 

uk  G  Vh,  Jk{uk)  =  inf  Jfc(u),  Jk(v)  =  Juk-i{v)  (21) 

vevh 

Since  Jk  is  not  Gateaux  differentiable,  one  can  use  (see  for  instance  [17,  18,  5])  the  function 
<Pe(x)  =  \J (x1  +  e2)  —  e  to  regularize  the  euclidian  norm  of  the  second  order  tensors.  With 
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this  technique,  the  material  is  not  completely  rigid  anymore,  and  so  it  is  difficult  to  capture 
accurately  the  shape  of  the  rigid  zone.  In  this  paper,  we  use  a  decomposition-coordination 
formulation  coupled  with  the  augmented  lagrangian  method  (see  [14,  13]).  This  approach  is 
accurate  in  detecting  the  visco-plastic  regions  and  permit  us  to  handle  the  locking  condition. 
Let 

Ah  =  Sh  G  L\V )3,  8h |t  G  Pi(T)3  VT  g  Th,  }, 

&h  =  { 0h ;  eh  g  L2(D),  eh\T  g  Pi(T)  VT  g  rA,  } 

and  the  Lagrangian  LI: 


Lk(v,cr,9,S)  =  /  r/|6|2  + 


/D 


fc— 1\ 


U 


k- 1 


ID 
■  Vt 


2X+-31le  + 


/  /e|<5|  -  p6  + 
Id  Jv 


u 


k- 1 12 


+  /  p(wfc_i  •  V)«fc_i  •  v  -  /  p/  •  u  + 


/x> 


/D 


(D'(v)  —  S)  :  a  +  {0  —  divu)  •  trace  cr 


(22) 


where  u  G  14,  d  G  A/,  represents  D'(v),  8  G  0^  stands  for  divu  and  the  lagrangian  multiplier 
cr  G  Ah  is  the  stress  inside  the  material.  Let  us  introduce  the  augmented  Lagrangian  Ck: 


£k(v,cr,9,8)  =  Lk(v,cr,9,8)  +  rD  f  \D\v)  -  5\2  +  rH  [  |0-divu|2,  (23) 

Jv  Jv 

where  rji  and  rv  are  strictly  positive  constants.  As  it  follows  from  [14,  13]  Ck  is  quadratic 
with  respect  with  v  and  its  saddle  points  coincide  with  those  of  Lk.  Thus,  the  minimization 
problem  (21)  becomes:  find  uk  G  14,  crk  G  A^,  8k  G  @h  and  Sk  G  Ah  such  that 

sup  Ck{uk ,  cr ,  9 , 8)  <  £k{uk ,  crk ,  6k ,  8k)  <  inf  Ck(v,crk  ,9k  ,8k).  (24) 

8.cr,S  v^Vh 


In  order  to  solve  the  above  saddle  point  problem  we  shall  use  an  Uzawa-type  algorithm 
(see  [14]).  For  this  let  us  put  <Tq_1  =  crk~l,9 q-1  =  9k~ 1  and  <5q-1  =  8k~1  and  let  us  introduce 
two  functions  fu  and  fo  useful  in  the  description  of  the  algorithm: 


fD(<r,a) 


1 

2 v(p)  +  a 


k{p) 

k'l 


a 


/ 


+ 


fH{cr,a) 


3  p(p)  +  tr  a 

< 

3A  (p)  +  2  r](p)  +  a 

0 

if  P  <  p* 
if  p  >  p* 


(25) 


(26) 


where  [a:]+  =  (x  +  \x\)/ 2  and  [x]_  =  (x  —  \x\)/ 2  represent  respectively  the  positive  and 
the  negative  part.  Let  us  remark  that,  by  taking  a  =  0,  (25)  and  (26)  give  us  the  complete 
model  of  our  material: 


D  =  fo(cr,  0)  + 


3 


Id  =  /(cr) 


At  the  iteration  i,  we  know  cr^_x , 8k_{,  9k_  1  and  we  compute  u f  ,cr 


k~  1  _ k— 1  rk—1  nk—1 


8k  1 , 9k  1  as  follows 


u^1  G  14, 


Ck{u~\cr:l,9ll,8:l)=  inf  £k(v,  0f_i  ,  Cl  )■ 

vevh 


<Sf_1  =  fD{<TiIi  +  2rDD(uk~1),2rD),  Ok~l  =  fH{crkZ^  +  2rHD(uk~1)),2rH) 

Kfc-1y  =  ((TiZly  +  ‘^rD(D'(uk-1)-Sk-1),  trace  erf”1  =  trace  <rkll  +  2rH{divuk~1  -0f_1). 

For  large  enough  i  =  imax  we  put  uk  =  ukzL,(rk  =  crkzL,Ok  =  OkzL  and  Sk  =  8kZ.L-  The 
interest  of  this  algorithm  is  that  it  transforms  the  non-differentiable  problem  into  a  sequence 
of  completely  standard  computations. 

We  shall  use  in  the  numerical  results  presented  in  the  next  sections  a  single  pass  (i.e. 
imax  _  2)  Uzawa-type  algorithm  to  deduce  the  following  numerical  approach  for  the  velocity 
problem. 

A) .  The  algorithm  starts  with  arbitraries  er°,  9 0  and  <5°. 

B) .  At  the  iteration  k.  we  have  <rk  1,  8k  1  and  9k  1 .  From  the  previous  density  problem 

we  know  also  the  density  and  therefore  the  distributions  of  p,r],  A  and  k.  Then  we  have 

to  compute  the  following  three  steps: 

1  —  First  Step: 

find  uk  €  Vh  such  that  Lk(uk:cr fc_1, <5fe_1)  =  inf  Ck(v,  crfc_1,  S^1).  (27) 

vevh 

This  problem  is  equivalent  to  find  uk  such  that: 

-dW(2rDD'(uk)  +  -rHdw(uk)I3)  =  div  {(Tk~l  -  2rD8k~ 1  -  - rH0k~1I3)  + 

3  3 

pfc-1(('Ufc-1  •  V),Mfe_1  —  /), 
with  appropriate  boundary  conditions. 

2  —  Second  Step:  Compute  explicitly  8k  and  8k  using  the  constitutive  equations  (25) 
and  (26): 

f  8k  =  fD(ak-1  +  2rDD(uk),2rD), 

1  9k  =  fH(*k-1  +  2rHD(uk),2rH). 

3  —  Third  Step:  Compute  crk  through  the  following  formulas: 

r  {aky  =  {(T^1)'  +  2rD(D'{vk)  -  8k), 

\  trace  ak  =  trace  erfc_1  +  2 vh  ( div  vk  —  . 


C).  The  algorithm  stops  when 


uk_uk  i||L2(p) 

Ill'll  VHP) 


is  small  enough. 


4.2  The  density  problem 

A  finite  volume  method  (see  for  instance  [12])  will  be  used  to  discretize  equation  (18).  Let 
denote  by  the  finite  volume  mesh,  which  is  given  by  a  family  of  disjoint  polygonal  connected 
subsets  of  M3  such  that  T>  is  the  union  of  the  closure  of  the  elements  of  h  is  the  greatest 
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diameter  of  a  control  volume  in  ICh-  The  finite  volume  mesh  ICh  is  obtained  from  a  finite 
element  triangulation  7/.  In  two  space  dimensions,  the  middle  of  each  side  of  a  triangle  T 
is  connected  to  the  center  of  mass  to  obtain  three  pieces.  A  control  volume  K  G  ICh,  which 
corresponds  to  the  node  A,  is  obtained  from  the  union  of  all  pieces  which  have  the  node  A. 

The  finite  volume  discrete  space  is  the  space  of  piecewise  constant  functions,  i.e.  we  are 
looking  for  the  solution  pfc  of  (18)  as  {p^k^  ^  £  /C h }•  In  order  to  simplify  the  notations  we 
shall  omit  in  all  this  subsection  the  indexes  h  and  k  and  pp  and  it/  will  be  denoted  simply  by 
p  and  u. 

Let  us  consider  two  control  volumes  K ,  P  G  ICh  with  a  common  interface  Ikl  =  Kf^P. 
Let  n  be  the  unit  normal  vector  to  Ikp  oriented  from  K  to  P.  The  we  define  the  flux  F(K,  P) 
at  the  interface  Ikp  as 

F(K,P)=  [  [u  ■  n]+. 

■J  Ikp 

Note  that  at  least  one  of  the  two  fluxes  F(K,P)  and  F(P,K )  is  vanishing.  If  we  denote  by 
A f{K)  the  set  of  all  neighbors  of  the  control  volume  K  then  the  finite  volume  numerical  scheme 
for  (18)  reads  as  a  linear  algebraic  system  for  the  unknowns  {pK)KeKh 

F(K,P)pk-F(P,K)pp  =  0,  for  all  K  G  Kh.  (28) 

PeM{K) 

If  a  volume  control  L  corresponds  to  a  node  which  is  on  the  boundary  dpV  then  we  put  pp  =  po 
and  we  eliminate  the  corresponding  equation. 


5  Penetration  in  concrete 


5.1  Material  parameters  model 


To  illustrate  the  predictive  capabilities  of  the  model,  its  application  to  concrete  is  presented. 
The  data  available  consists  of  laboratory  quasi-static  unconfined  and  confined  compression 
tests  for  confining  pressures  in  the  range  50-450  MPa  under  a  strain  rate  of  10-6/s  and  both 
confined  and  unconfined  Split-Hopkinson  bar  data  at  strain  rates  of  60/s  to  160/s.  The  ambient 
density  of  the  material  is:po  =  2000kgm“3  and  the  unconfined  strength  is  of  60.00  MPa(data 
after  Schmdit[29].  The  locking  density  is  of  p*  =  2,600  kgm-3  and  locking  pressure  p*  =  0.5 
GPa. 

The  pressure  -  density  relationship  for  this  material  can  be  approximated  by 


p(p) 

*  P  / 1  Po  \ 

=  p  .  1  b 

Ki 

TS 

o 

IA 

A 

* 

T5 

1 

TS 

o 

p{p) 

>P*, 

if  p>  p* 

(29) 


The  following  law  of  variation  of  the  yield  limit  with  the  current  density  describes  well  the 
data: 

k(p)  =  k0  +  P(1  -  — ) 

P 

where  /3  =  (k*  —  kq)po/{p*  ~  Po)  and  Ko,  are  the  yield  stress  corresponding  to  the  density 
of  the  undeformed  (p  =  po)  and  locked  medium  (p  =  p*).  In  the  numerical  simulations,  we 
have  taken  kq  =  100 MPa,  and  k*  =  800 MPa. 
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For  the  viscosity  coefficients  r)(p)  and  A (p)  a  similar  variation  law  as  for  the  yield  limit  is 
used 

v{p)  =  Vo  +  7(1  -  —  f,  Kp)  =  Ao  +  <5(1  -  —  f, 

P  P 

where  7  =  (rj*  —  Vo)(po/(p*  ~  Po ))2  and  <5  =  (A*  —  Ao)(po/(p*  —  /°o))2-  The  set  of  values  chosen 
in  the  numerical  simulations  are:  70  =  20.  kPa-s, 77*  =  5.  kPa-s,  Ao  =  10  kPa-s  and  A*  =  1.0 
MPa-s. 


5.2  Boundary  conditions 

On  the  boundary  do T>,  which  consists  of  the  polygonal  line  OFEC,  the  velocity  is  Vez,  where 
V  is  the  impact  velocity.  Simulations  were  performed  for  an  impact  velocity  of  500  ms  1 .  Since 
the  control  volume  T>  does  not  touch  the  tunnel  we  have  d\D  =  0.  Preliminary  simulations,  in 
which  a  part  of  the  tunnel  was  modelled  have  shown  that  the  material  in  V  which  is  behind 
the  projectile  is  rigid  and  thus  does  not  influence  the  computations. 

The  frictional  contact  boundary  is  the  boundary  of  the  projectile,  i.e.  ABC  on  the  figure 
1.  We  have  chosen  Lim-Ashby-Klepaczko  model  (see  [20])  to  describe  the  influence  of  the 
velocity  on  the  friction  coefficient  between  the  concrete  target  and  the  metallic  penetrator 

p{v)  =  a[l  -  b  log ( 1  +  — )] 

where  a  =  0.5,  b  =  0.263  and  vq  =  1.0  ms-1.  The  boundary  condition  un  =  0  prescribed  on  the 
projectile  is  accurate  everywhere  apart  from  two  small  zones.  The  first  one  is  located  on  the 
nose  of  the  projectile  in  the  neighborhood  of  A.  In  this  zone,  rupture  in  the  target  material  and 
creation  of  a  free  boundary,  which  begins  somewhere  on  OA  and  ends  somewhere  on  AB ,  but 
still  very  close  to  A  is  to  be  expected.  In  the  second  zone,  which  is  located  on  the  body  of  the 
projectile  behind  and  near  B,  the  target  material  is  no  longer  in  contact  with  the  projectile. 
This  corresponds  to  a  second  free  boundary  zone  which  starts  and  ends  somewhere  on  the 
body  BC.  The  formulation  and  solution  of  these  two  free  boundary  problems  is  beyond  the 
scope  of  this  study.  However,  the  influence  of  these  two  zones  on  the  resistance  to  penetration 
seems  to  be  very  small. 

5.3  Mesh  adaptation 

In  order  to  capture  sharply  the  shape  of  the  visco-plactic  zone,  we  have  used  the  anisotropic 
mesh  generator  BAMG  (see  [15,  16]).  This  generator  requires  a  governing  field  for  re-meshing 
and  refines  the  zones  with  high  second  derivative.  The  field  E  we  used  is  the  square  root  of 
the  dissipative  energy  associated  to  a  velocity  field  u  computed  through  <5  and  9  defined  in 
(22).  _ 

E=  y/tr(u)  :  D(u)  =  \J r](p)\6\2  +  k(P)\S\  -  p(p)9  +  +  3A(d)  Q  (30) 

Note  that  E  is  continuous,  but  its  first  derivative  is  discontinuous  across  the  zones  of  contact 
between  the  rigid  and  visco-plastic  regions;  that  means  that  the  second  derivative  is  very  large 
on  this  surface  that  the  mesh  generator  will  emphasize  this  boundary.  Other  regions  have  to 
be  refined,  those  ones  with  high  rates  of  deformation.  As  we  can  expect,  the  zone  around  the 
nose  of  the  projectile  will  be  also  privileged. 

In  figure  2  we  have  plotted  the  initial  finite  element  mesh  of  the  domain  OABCEF,  and  the 
distribution  of  the  deformation  rate  \D(u)\.  Note  that  the  boundary  between  the  rigid  zone 
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and  the  visco-plastic  one  is  not  well  captured.  This  mesh  has  2304  nodes  and  4356  triangular 
elements.  The  finite  element  mesh  obtained  after  remeshing  is  plotted  in  Figure  3.  Using  6992 
nodes  and  13,719  triangular  elements  this  mesh  is  able  to  give  a  satisfactory  description  of  the 
rigid  zone. 


Z-axis 


Z-axis 


Figure  2:  The  initial  finite  element  mesh  (left),  and  the  distribution  of  the  deformation  rate 
|Z)(w)|  (right).  Note  that  the  boundary  between  the  rigid  zone  and  the  visco-plastic  one  is  not 
well  captured. 


Figure  3:  The  finite  element  mesh  obtained  after  remeshing  (left),  and  the  distribution  of 
the  deformation  rate  \D(u)\  (right).  Note  that  the  boundary  between  the  rigid  zone  and  the 
visco-plastic  one  is  sharply  captured. 
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6  Results  and  discussion 


Figure  4  shows  the  distribution  of  density  in  the  target.  Note  the  existence  of  three  distinct 
zones:  (1)  a  fully  compacted  zone  around  the  penetrator,  where  the  density  is  everywhere 
equal  to  the  locking  density  p*  ,  (2)  a  compacted  where  the  density  varies  between  po,  the 
density  of  the  intact  material,  and  p*;  (3)  a  zone  that  remains  unaffected  by  the  impact  event 
where  the  density  is  everywhere  equal  to  po-  The  streamlines  are  shown  in  figure  streamlines. 
It  is  worth  noticing  that  all  the  particles  ahead  of  the  projectile  which  are  at  a  distance  from 
the  centerline  less  than  R  will  enter  the  fully  compacted  zone  (1). 


Z-axis 


Figure  4:  The  distribution  of  density  in  the  target  and  the  streamlines  computed  in  Kg  m  . 

In  figure  5  it  is  shown  the  distribution  of  the  second  invariant  of  the  rate  of  deformation 
deviator  \D\u)\.  It  is  seen  that  a  zone  of  intense  plastic  deformation  develops  around  the 
penetrator  and  it  extends  outward  to  approximatively  3  projectile  radii  from  the  centerline. 
The  maximum  rate  of  deformation  is  achieved  on  the  nose  tip.  Outside  this  zone,  | D'(u)\ 
vanishes,  hence  according  to  the  model  no  plastic  deformation  can  occur. 
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Z-axis 


Figure  5:  The  distribution  of  the  second  invariant  of  the  rate  of  deformation  deviator  \D'(u)\ 
computed  in  s_1. 

The  distribution  of  the  rate  of  volumetric  deformation  divu  is  displayed  in  Figure  6.  We 
note  that  around  the  penetrator,  a  fully  compacted  state  is  achieved,  the  maximum  com¬ 
paction  being  in  the  nose  zone.  Because  the  prescribed  boundary  conditions  do  not  reflect 
that  the  target  material  is  not  in  contact  with  the  projectile  on  a  limited  zone  of  the  body 
of  the  projectile,  divit  is  slightly  positive  in  a  small  region  behind  the  projectile.  This  zone 
corresponds  to  positive  normal  stresses.  This  can  be  seen  in  Figure  7  where  it  is  shown  the 
distribution  of  the  normal  stress  crn  along  the  projectile.  High  compressive  values  on  the  tip  of 
the  nose  and  very  small  values  on  the  shank  are  obtained.  The  tangential  velocity  ut  on  the 
projectile,  plotted  in  Figure  7,  is  very  small  on  the  tip  and  gradually  increases  until  it  reaches 
the  value  of  the  impact  velocity  on  the  shank. 
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Figure  6:  The  distribution  of  the  rate  of  volumetric  deformation  diva  computed  in  s 


Figure  7:  Left:  the  distribution  of  tangential  velocity  ut  along  the  projectile  computed  in 
ms-1.  Right  :  The  distribution  of  the  normal  stress  crn  along  the  projectile  computed  in  Pa. 

Let  us  analyze  now  the  stress  distribution  in  the  target  and  determine  the  zones  where 
tensile  failure  may  occur.  It  is  assumed  that  concrete  tensile  failure  can  be  modelled  with  the 
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classic  maximum  tensile  strength  criterion  (see  [25]),  i.e.  the  material  fails  if  the  maximum 
of  the  principal  stresses  reaches  a  critical  limit  denoted  /.  The  tensile  strength  depends  on 
the  level  of  compaction  of  the  material.  A  power  law  variation  of  this  tensile  strength  of  the 
material  with  the  current  density  is  assumed 

f{p)  =  fo  T  o(l  -  — )2 
P 

where  a  =  (/*  —  fo)(po/(p*  —  Po))2  and  fo,  f*  are  the  tensile  strength  corresponding  to  the 
density  of  the  non  deformed  (p  =  po)  and  locked  medium  (p  =  p*).  For  the  numerical  tests 
we  have  used  fo  =  33.3 MPa,  and  f*  =  233.3 MPa.  The  plot  of  the  distribution  of  agg/f  is 
plotted  in  figure  8.  We  remark  that  in  a  wide  region  ahead  of  the  projectile  the  tensile  failure 
can  occur  along  Ozr  planes  orthogonal  to  egg  (see  figure  1).  Following  the  stream  lines  a  part 
of  this  region  will  be  fully  compacted,  hence  in  the  lagrangian  configuration  fracture  will  not 
be  observed  near  the  projectile. 

From  the  analysis  of  the  spatial  distribution  of  damage  in  the  target  one  can  gain  funda¬ 
mental  understanding  of  the  stability  of  the  penetrator  trajectory.  We  have  found  that  damage 
occurs  ahead  the  projectile  and  along  planes  which  can  be  not  symmetric  with  respect  to  the 
penetrator  centerline.  This  damage  induced  anisotropy  can  explain  trajectory  deviations. 
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Figure  8:  The  distribution  of  crgg/f  in  the  target. 
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7  Conclusions 


A  model  for  describing  steady-state  penetration  of  a  rigid  penetrator  into  a  geologic  or  geo¬ 
logically  derived  material  was  proposed. 

Since  for  low-to-intermediate  impact  velocities,  the  impacted  medium  displays  both  solid¬ 
like  and  fluid-like  behavior  a  rigid  visco-plastic  fluid  type  constitutive  was  developed.  To 
capture  the  observed  dependence  of  yielding  and  strain-rate  sensitivity  on  compaction,  an  ex¬ 
plicit  dependence  of  the  yield  limit  and  viscosity  coefficients  on  density  was  considered.  To 
reflect  the  observed  characteristics  of  the  pressure- density  relationships  in  geologic  materials 
when  subjected  to  high-pressure  (of  the  order  of  GPa)  the  hypothesis  of  rigid  unloading  and 
locking  medium  was  adopted  (i.e.  the  density  cannot  exceed  a  critical  limit).  Variational  for¬ 
mulations  and  algorithms  for  solving  the  minimization  problems  for  the  velocity  field  using  the 
finite  element  method  were  developed  while  finite-volume  techniques  were  adopted  for  solving 
the  hyperbolic  mass  conservation  equation.  To  solve  the  velocity  problem  a  decomposition- 
coordination  formulation  coupled  with  the  augmented  lagrangian  method  is  used.  The  model 
was  applied  to  penetration  into  concrete.  The  deformation  behavior  was  described  using  the 
viscoplastic  compressible  fluid  model  developed  while  fracture  was  modeled  using  the  classic 
maximum  tensile  stress  criterion.  The  model  predicts  that  around  the  penetrator,  a  fully  com¬ 
pacted  state  is  achieved,  the  maximum  compaction  being  in  the  nose  zone.  Fracture  occurs 
ahead  of  the  penetrator  and  along  planes  which  can  be  not  symmetric  with  respect  to  the 
penetrator  centerline.  This  induced  anisotropy  can  explain  trajectory  deviations. 
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